home *** CD-ROM | disk | FTP | other *** search
/ The Atari Compendium / The Atari Compendium (Toad Computers) (1994).iso / files / prgtools / gnustuff / tos / updates / update36.zoo / libg++ / iosrc / floatconv.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-11-28  |  71.0 KB  |  2,341 lines
  1. /* 
  2. Copyright (C) 1993 Free Software Foundation
  3.  
  4. This file is part of the GNU IO Library.  This library is free
  5. software; you can redistribute it and/or modify it under the
  6. terms of the GNU General Public License as published by the
  7. Free Software Foundation; either version 2, or (at your option)
  8. any later version.
  9.  
  10. This library is distributed in the hope that it will be useful,
  11. but WITHOUT ANY WARRANTY; without even the implied warranty of
  12. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13. GNU General Public License for more details.
  14.  
  15. You should have received a copy of the GNU General Public License
  16. along with GNU CC; see the file COPYING.  If not, write to
  17. the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  
  19. As a special exception, if you link this library with files
  20. compiled with a GNU compiler to produce an executable, this does not cause
  21. the resulting executable to be covered by the GNU General Public License.
  22. This exception does not however invalidate any other reasons why
  23. the executable file might be covered by the GNU General Public License. */
  24.  
  25. #include <libioP.h>
  26. #ifdef USE_DTOA
  27. /****************************************************************
  28.  *
  29.  * The author of this software is David M. Gay.
  30.  *
  31.  * Copyright (c) 1991 by AT&T.
  32.  *
  33.  * Permission to use, copy, modify, and distribute this software for any
  34.  * purpose without fee is hereby granted, provided that this entire notice
  35.  * is included in all copies of any software which is or includes a copy
  36.  * or modification of this software and in all copies of the supporting
  37.  * documentation for such software.
  38.  *
  39.  * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
  40.  * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
  41.  * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
  42.  * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
  43.  *
  44.  ***************************************************************/
  45.  
  46. /* Some cleaning up by Per Bothner, bothner@cygnus.com, 1992, 1993.
  47.    Re-written to not need static variables
  48.    (except result, result_k, HIWORD, LOWORD). */
  49.  
  50. /* Please send bug reports to
  51.         David M. Gay
  52.         AT&T Bell Laboratories, Room 2C-463
  53.         600 Mountain Avenue
  54.         Murray Hill, NJ 07974-2070
  55.         U.S.A.
  56.         dmg@research.att.com or research!dmg
  57.  */
  58.  
  59. /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
  60.  *
  61.  * This strtod returns a nearest machine number to the input decimal
  62.  * string (or sets errno to ERANGE).  With IEEE arithmetic, ties are
  63.  * broken by the IEEE round-even rule.  Otherwise ties are broken by
  64.  * biased rounding (add half and chop).
  65.  *
  66.  * Inspired loosely by William D. Clinger's paper "How to Read Floating
  67.  * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
  68.  *
  69.  * Modifications:
  70.  *
  71.  *      1. We only require IEEE, IBM, or VAX double-precision
  72.  *              arithmetic (not IEEE double-extended).
  73.  *      2. We get by with floating-point arithmetic in a case that
  74.  *              Clinger missed -- when we're computing d * 10^n
  75.  *              for a small integer d and the integer n is not too
  76.  *              much larger than 22 (the maximum integer k for which
  77.  *              we can represent 10^k exactly), we may be able to
  78.  *              compute (d*10^k) * 10^(e-k) with just one roundoff.
  79.  *      3. Rather than a bit-at-a-time adjustment of the binary
  80.  *              result in the hard case, we use floating-point
  81.  *              arithmetic to determine the adjustment to within
  82.  *              one bit; only in really hard cases do we need to
  83.  *              compute a second residual.
  84.  *      4. Because of 3., we don't need a large table of powers of 10
  85.  *              for ten-to-e (just some small tables, e.g. of 10^k
  86.  *              for 0 <= k <= 22).
  87.  */
  88.  
  89. /*
  90.  * #define IEEE_8087 for IEEE-arithmetic machines where the least
  91.  *      significant byte has the lowest address.
  92.  * #define IEEE_MC68k for IEEE-arithmetic machines where the most
  93.  *      significant byte has the lowest address.
  94.  * #define Sudden_Underflow for IEEE-format machines without gradual
  95.  *      underflow (i.e., that flush to zero on underflow).
  96.  * #define IBM for IBM mainframe-style floating-point arithmetic.
  97.  * #define VAX for VAX-style floating-point arithmetic.
  98.  * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
  99.  * #define No_leftright to omit left-right logic in fast floating-point
  100.  *      computation of dtoa.
  101.  * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
  102.  * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
  103.  *      that use extended-precision instructions to compute rounded
  104.  *      products and quotients) with IBM.
  105.  * #define ROUND_BIASED for IEEE-format with biased rounding.
  106.  * #define Inaccurate_Divide for IEEE-format with correctly rounded
  107.  *      products but inaccurate quotients, e.g., for Intel i860.
  108.  * #define KR_headers for old-style C function headers.
  109.  */
  110.  
  111. #ifdef DEBUG
  112. #include <stdio.h>
  113. #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
  114. #endif
  115.  
  116. #ifdef __STDC__
  117. #include <stdlib.h>
  118. #include <string.h>
  119. #include <float.h>
  120. #define CONST const
  121. #else
  122. #define CONST
  123. #define KR_headers
  124.  
  125. /* In this case, we assume IEEE floats. */
  126. #define FLT_ROUNDS 1
  127. #define FLT_RADIX 2
  128. #define DBL_DIG 15
  129. #define DBL_MAX_10_EXP 308
  130. #define DBL_MAX_EXP 1024
  131. #endif
  132.  
  133. #include <errno.h>
  134. #ifndef __MATH_H__
  135. #include <math.h>
  136. #endif
  137.  
  138. #ifdef Unsigned_Shifts
  139. #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000;
  140. #else
  141. #define Sign_Extend(a,b) /*no-op*/
  142. #endif
  143.  
  144. #if defined(__i386__) || defined(clipper) || defined(MIPSEL)
  145. #define IEEE_8087
  146. #endif
  147.  
  148. #if defined(__sparc__) || defined(sparc) || defined(MIPSEB) || defined(atarist)
  149. #define IEEE_MC68k
  150. #endif
  151.  
  152. #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
  153.  
  154. #if FLT_RADIX==16
  155. #define IBM
  156. #else
  157. #if DBL_MANT_DIG==56
  158. #define VAX
  159. #else
  160. #if DBL_MANT_DIG==53 && DBL_MAX_10_EXP==308
  161. #define IEEE_Unknown
  162. #else
  163. Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
  164. #endif
  165. #endif
  166. #endif
  167. #endif
  168.  
  169. #ifdef IEEE_8087
  170. #define HIWORD 1
  171. #define LOWORD 0
  172. #define TEST_ENDIANNESS  /* nothing */
  173. #else
  174. #if defined(IEEE_MC68k)
  175. #define HIWORD 0
  176. #define LOWORD 1
  177. #define TEST_ENDIANNESS  /* nothing */
  178. #else
  179. static int HIWORD = -1, LOWORD;
  180. static void test_endianness()
  181. {
  182.     union doubleword {
  183.     double d;
  184.     unsigned long u[2];
  185.     } dw;
  186.     dw.d = 10;
  187.     if (dw.u[0] != 0) /* big-endian */
  188.     HIWORD=0, LOWORD=1;
  189.     else
  190.     HIWORD=1, LOWORD=0;
  191. }
  192. #define TEST_ENDIANNESS  if (HIWORD<0) test_endianness();
  193. #endif
  194. #endif
  195.  
  196. #if 0
  197. union {
  198.   double d;
  199.   unsigned long x[2];
  200. } _temp;
  201. #endif
  202. #define word0(x) ((unsigned long *)&x)[HIWORD]
  203. #if 0
  204. #define word0(X) (_temp.d = X, _temp.x[HIWORD])
  205. #define setword0(D,X) (_temp.d = D, _temp.x[HIWORD] = X, D = _temp.d)
  206. #endif
  207. #define word1(x) ((unsigned long *)&x)[LOWORD]
  208.  
  209. /* The following definition of Storeinc is appropriate for MIPS processors. */
  210. #if defined(IEEE_8087) + defined(VAX)
  211. #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
  212. ((unsigned short *)a)[0] = (unsigned short)c, a++)
  213. #else
  214. #if defined(IEEE_MC68k)
  215. #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
  216. ((unsigned short *)a)[1] = (unsigned short)c, a++)
  217. #else
  218. #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
  219. #endif
  220. #endif
  221.  
  222. /* #define P DBL_MANT_DIG */
  223. /* Ten_pmax = floor(P*log(2)/log(5)) */
  224. /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
  225. /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
  226. /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
  227.  
  228. #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(IEEE_Unknown)
  229. #define Exp_shift  20
  230. #define Exp_shift1 20
  231. #define Exp_msk1    0x100000
  232. #define Exp_msk11   0x100000
  233. #define Exp_mask  0x7ff00000
  234. #define P 53
  235. #define Bias 1023
  236. #define IEEE_Arith
  237. #define Emin (-1022)
  238. #define Exp_1  0x3ff00000
  239. #define Exp_11 0x3ff00000
  240. #define Ebits 11
  241. #define Frac_mask  0xfffff
  242. #define Frac_mask1 0xfffff
  243. #define Ten_pmax 22
  244. #define Bletch 0x10
  245. #define Bndry_mask  0xfffff
  246. #define Bndry_mask1 0xfffff
  247. #define LSB 1
  248. #define Sign_bit 0x80000000
  249. #define Log2P 1
  250. #define Tiny0 0
  251. #define Tiny1 1
  252. #define Quick_max 14
  253. #define Int_max 14
  254. #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */
  255. #else
  256. #undef  Sudden_Underflow
  257. #define Sudden_Underflow
  258. #ifdef IBM
  259. #define Exp_shift  24
  260. #define Exp_shift1 24
  261. #define Exp_msk1   0x1000000
  262. #define Exp_msk11  0x1000000
  263. #define Exp_mask  0x7f000000
  264. #define P 14
  265. #define Bias 65
  266. #define Exp_1  0x41000000
  267. #define Exp_11 0x41000000
  268. #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
  269. #define Frac_mask  0xffffff
  270. #define Frac_mask1 0xffffff
  271. #define Bletch 4
  272. #define Ten_pmax 22
  273. #define Bndry_mask  0xefffff
  274. #define Bndry_mask1 0xffffff
  275. #define LSB 1
  276. #define Sign_bit 0x80000000
  277. #define Log2P 4
  278. #define Tiny0 0x100000
  279. #define Tiny1 0
  280. #define Quick_max 14
  281. #define Int_max 15
  282. #else /* VAX */
  283. #define Exp_shift  23
  284. #define Exp_shift1 7
  285. #define Exp_msk1    0x80
  286. #define Exp_msk11   0x800000
  287. #define Exp_mask  0x7f80
  288. #define P 56
  289. #define Bias 129
  290. #define Exp_1  0x40800000
  291. #define Exp_11 0x4080
  292. #define Ebits 8
  293. #define Frac_mask  0x7fffff
  294. #define Frac_mask1 0xffff007f
  295. #define Ten_pmax 24
  296. #define Bletch 2
  297. #define Bndry_mask  0xffff007f
  298. #define Bndry_mask1 0xffff007f
  299. #define LSB 0x10000
  300. #define Sign_bit 0x8000
  301. #define Log2P 1
  302. #define Tiny0 0x80
  303. #define Tiny1 0
  304. #define Quick_max 15
  305. #define Int_max 15
  306. #endif
  307. #endif
  308.  
  309. #ifndef IEEE_Arith
  310. #define ROUND_BIASED
  311. #endif
  312.  
  313. #ifdef RND_PRODQUOT
  314. #define rounded_product(a,b) a = rnd_prod(a, b)
  315. #define rounded_quotient(a,b) a = rnd_quot(a, b)
  316. extern double rnd_prod(double, double), rnd_quot(double, double);
  317. #else
  318. #define rounded_product(a,b) a *= b
  319. #define rounded_quotient(a,b) a /= b
  320. #endif
  321.  
  322. #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
  323. #define Big1 0xffffffff
  324.  
  325. #define Kmax 15
  326.  
  327. /* (1<<BIGINT_MINIMUM_K) is the minimum number of words to allocate
  328.    in a Bigint.  dtoa usually manages with 1<<2, and has not been
  329.    known to need more than 1<<3.  */
  330.  
  331. #define BIGINT_MINIMUM_K 3
  332.  
  333. struct Bigint {
  334.   struct Bigint *next;
  335.   int k;        /* Parameter given to Balloc(k) */
  336.   int maxwds;        /* Allocated space: equals 1<<k. */
  337.   short on_stack;    /* 1 if stack-allocated. */
  338.   short sign;        /* 0 if value is positive or zero; 1 if negative. */
  339.   int wds;        /* Current length. */
  340.   unsigned long x[1<<BIGINT_MINIMUM_K]; /* Actually: x[maxwds] */
  341. };
  342.  
  343. #define BIGINT_HEADER_SIZE \
  344.   (sizeof(Bigint) - (1<<BIGINT_MINIMUM_K) * sizeof(unsigned long))
  345.  
  346. typedef struct Bigint Bigint;
  347.  
  348. /* Initialize a stack-allocated Bigint. */
  349.  
  350. Bigint *
  351. Binit
  352. #ifdef KR_headers
  353.         (v) Bigint *v;
  354. #else
  355.         (Bigint *v)
  356. #endif
  357. {
  358.   v->on_stack = 1;
  359.   v->k = BIGINT_MINIMUM_K;
  360.   v->maxwds = 1 << BIGINT_MINIMUM_K;
  361.   v->sign = v->wds = 0;
  362.   return v;
  363. }
  364.  
  365. /* Allocate a Bigint with '1<<k' big digits. */
  366.  
  367. static Bigint *
  368. Balloc
  369. #ifdef KR_headers
  370.         (k) int k;
  371. #else
  372.         (int k)
  373. #endif
  374. {
  375.   int x;
  376.   Bigint *rv;
  377.  
  378.   if (k < BIGINT_MINIMUM_K)
  379.     k = BIGINT_MINIMUM_K;
  380.  
  381.   x = 1 << k;
  382.   rv = (Bigint *)
  383.     malloc(BIGINT_HEADER_SIZE + x * sizeof(unsigned long));
  384.   rv->k = k;
  385.   rv->maxwds = x;
  386.   rv->sign = rv->wds = 0;
  387.   rv->on_stack = 0;
  388.   return rv;
  389. }
  390.  
  391. static void
  392. Bfree
  393. #ifdef KR_headers
  394.         (v) Bigint *v;
  395. #else
  396.         (Bigint *v)
  397. #endif
  398. {
  399.   if (v && !v->on_stack)
  400.     free (v);
  401. }
  402.  
  403. static void
  404. Bcopy
  405. #ifdef KR_headers
  406.         (x, y) Bigint *x, *y;
  407. #else
  408.         (Bigint *x, Bigint *y)
  409. #endif
  410. {
  411.   register unsigned long *xp, *yp;
  412.   register int i = y->wds;
  413.   x->sign = y->sign;
  414.   x->wds = i;
  415.   for (xp = x->x, yp = y->x; --i >= 0; )
  416.     *xp++ = *yp++;
  417. }
  418.  
  419. /* Make sure b has room for at least 1<<k big digits. */
  420.  
  421. static Bigint *
  422. Brealloc
  423. #ifdef KR_headers
  424.         (b, k) Bigint *b; int k;
  425. #else
  426.         (Bigint * b, int k)
  427. #endif
  428. {
  429.   if (b == NULL)
  430.     return Balloc(k);
  431.   if (b->k >= k)
  432.     return b;
  433.   else
  434.     {
  435.       Bigint *rv = Balloc (k);
  436.       Bcopy(rv, b);
  437.       Bfree(b);
  438.       return rv;
  439.     }
  440. }
  441.  
  442. /* Return b*m+a.  b is modified.
  443.    Assumption:  0xFFFF*m+a fits in 32 bits. */
  444.  
  445. static Bigint *
  446. multadd
  447. #ifdef KR_headers
  448.         (b, m, a) Bigint *b; int m, a;
  449. #else
  450.         (Bigint *b, int m, int a)
  451. #endif
  452. {
  453.         int i, wds;
  454.         unsigned long *x, y;
  455.         unsigned long xi, z;
  456.  
  457.         wds = b->wds;
  458.         x = b->x;
  459.         i = 0;
  460.         do {
  461.                 xi = *x;
  462.                 y = (xi & 0xffff) * m + a;
  463.                 z = (xi >> 16) * m + (y >> 16);
  464.                 a = (int)(z >> 16);
  465.                 *x++ = (z << 16) + (y & 0xffff);
  466.                 }
  467.                 while(++i < wds);
  468.         if (a) {
  469.                 if (wds >= b->maxwds)
  470.                         b = Brealloc(b, b->k+1);
  471.                 b->x[wds++] = a;
  472.                 b->wds = wds;
  473.                 }
  474.         return b;
  475.         }
  476.  
  477. static Bigint *
  478. s2b
  479. #ifdef KR_headers
  480.         (result, s, nd0, nd, y9)
  481.     Bigint *result; CONST char *s; int nd0, nd; unsigned long y9;
  482. #else
  483.         (Bigint *result, CONST char *s, int nd0, int nd, unsigned long y9)
  484. #endif
  485. {
  486.   int i, k;
  487.   long x, y;
  488.  
  489.   x = (nd + 8) / 9;
  490.   for(k = 0, y = 1; x > y; y <<= 1, k++) ;
  491.   result = Brealloc(result, k);
  492.   result->x[0] = y9;
  493.   result->wds = 1;
  494.  
  495.   i = 9;
  496.   if (9 < nd0)
  497.     {
  498.       s += 9;
  499.       do
  500.     result = multadd(result, 10, *s++ - '0');
  501.       while (++i < nd0);
  502.       s++;
  503.     }
  504.   else
  505.     s += 10;
  506.   for(; i < nd; i++)
  507.     result = multadd(result, 10, *s++ - '0');
  508.   return result;
  509. }
  510.  
  511. static int
  512. hi0bits
  513. #ifdef KR_headers
  514.         (x) register unsigned long x;
  515. #else
  516.         (register unsigned long x)
  517. #endif
  518. {
  519.         register int k = 0;
  520.  
  521.         if (!(x & 0xffff0000)) {
  522.                 k = 16;
  523.                 x <<= 16;
  524.                 }
  525.         if (!(x & 0xff000000)) {
  526.                 k += 8;
  527.                 x <<= 8;
  528.                 }
  529.         if (!(x & 0xf0000000)) {
  530.                 k += 4;
  531.                 x <<= 4;
  532.                 }
  533.         if (!(x & 0xc0000000)) {
  534.                 k += 2;
  535.                 x <<= 2;
  536.                 }
  537.         if (!(x & 0x80000000)) {
  538.                 k++;
  539.                 if (!(x & 0x40000000))
  540.                         return 32;
  541.                 }
  542.         return k;
  543.         }
  544.  
  545. static int
  546. lo0bits
  547. #ifdef KR_headers
  548.         (y) unsigned long *y;
  549. #else
  550.         (unsigned long *y)
  551. #endif
  552. {
  553.         register int k;
  554.         register unsigned long x = *y;
  555.  
  556.         if (x & 7) {
  557.                 if (x & 1)
  558.                         return 0;
  559.                 if (x & 2) {
  560.                         *y = x >> 1;
  561.                         return 1;
  562.                         }
  563.                 *y = x >> 2;
  564.                 return 2;
  565.                 }
  566.         k = 0;
  567.         if (!(x & 0xffff)) {
  568.                 k = 16;
  569.                 x >>= 16;
  570.                 }
  571.         if (!(x & 0xff)) {
  572.                 k += 8;
  573.                 x >>= 8;
  574.                 }
  575.         if (!(x & 0xf)) {
  576.                 k += 4;
  577.                 x >>= 4;
  578.                 }
  579.         if (!(x & 0x3)) {
  580.                 k += 2;
  581.                 x >>= 2;
  582.                 }
  583.         if (!(x & 1)) {
  584.                 k++;
  585.                 x >>= 1;
  586.                 if (!x & 1)
  587.                         return 32;
  588.                 }
  589.         *y = x;
  590.         return k;
  591.         }
  592.  
  593. static Bigint *
  594. i2b
  595. #ifdef KR_headers
  596.         (result, i) Bigint *result; int i;
  597. #else
  598.         (Bigint* result, int i)
  599. #endif
  600. {
  601.   result = Brealloc(result, 1);
  602.   result->x[0] = i;
  603.   result->wds = 1;
  604.   return result;
  605. }
  606.  
  607. /* Do: c = a * b. */
  608.  
  609. static Bigint *
  610. mult
  611. #ifdef KR_headers
  612.         (c, a, b) Bigint *a, *b, *c;
  613. #else
  614.         (Bigint *c, Bigint *a, Bigint *b)
  615. #endif
  616. {
  617.         int k, wa, wb, wc;
  618.         unsigned long carry, y, z;
  619.         unsigned long *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
  620.         unsigned long z2;
  621.         if (a->wds < b->wds) {
  622.                 Bigint *tmp = a;
  623.                 a = b;
  624.                 b = tmp;
  625.                 }
  626.         k = a->k;
  627.         wa = a->wds;
  628.         wb = b->wds;
  629.         wc = wa + wb;
  630.         if (wc > a->maxwds)
  631.                 k++;
  632.     c = Brealloc(c, k);
  633.         for(x = c->x, xa = x + wc; x < xa; x++)
  634.                 *x = 0;
  635.         xa = a->x;
  636.         xae = xa + wa;
  637.         xb = b->x;
  638.         xbe = xb + wb;
  639.         xc0 = c->x;
  640.         for(; xb < xbe; xb++, xc0++) {
  641.                 if (y = *xb & 0xffff) {
  642.                         x = xa;
  643.                         xc = xc0;
  644.                         carry = 0;
  645.                         do {
  646.                                 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
  647.                                 carry = z >> 16;
  648.                                 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
  649.                                 carry = z2 >> 16;
  650.                                 Storeinc(xc, z2, z);
  651.                                 }
  652.                                 while(x < xae);
  653.                         *xc = carry;
  654.                         }
  655.                 if (y = *xb >> 16) {
  656.                         x = xa;
  657.                         xc = xc0;
  658.                         carry = 0;
  659.                         z2 = *xc;
  660.                         do {
  661.                                 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
  662.                                 carry = z >> 16;
  663.                                 Storeinc(xc, z, z2);
  664.                                 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
  665.                                 carry = z2 >> 16;
  666.                                 }
  667.                                 while(x < xae);
  668.                         *xc = z2;
  669.                         }
  670.                 }
  671.         for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
  672.         c->wds = wc;
  673.         return c;
  674.         }
  675.  
  676. /* Returns b*(5**k).  b is modified. */
  677. /* Re-written by Per Bothner to not need a static list. */
  678.  
  679. static Bigint *
  680. pow5mult
  681. #ifdef KR_headers
  682.         (b, k) Bigint *b; int k;
  683. #else
  684.         (Bigint *b, int k)
  685. #endif
  686. {
  687.   static int p05[6] = { 5, 25, 125, 625, 3125, 15625 };
  688.  
  689.   for (; k > 6; k -= 6)
  690.     b = multadd(b, 15625, 0); /* b *= 5**6 */
  691.   if (k == 0)
  692.     return b;
  693.   else
  694.     return multadd(b, p05[k-1], 0);
  695. }
  696.  
  697. /* Re-written by Per Bothner so shift can be in place. */
  698.  
  699. static Bigint *
  700. lshift
  701. #ifdef KR_headers
  702.     (b, k) Bigint *b; int k;
  703. #else
  704.         (Bigint *b, int k)
  705. #endif
  706. {
  707.   int i;
  708.   unsigned long *x, *x1, *xe;
  709.   int old_wds = b->wds;
  710.   int n = k >> 5;
  711.   int k1 = b->k;
  712.   int n1 = n + old_wds + 1;
  713.  
  714.   if (k == 0)
  715.     return b;
  716.  
  717.   for(i = b->maxwds; n1 > i; i <<= 1)
  718.     k1++;
  719.   b = Brealloc(b, k1);
  720.  
  721.   xe = b->x; /* Source limit */
  722.   x = xe + old_wds; /* Source pointer */
  723.   x1 = x + n; /* Destination pointer */
  724.   if (k &= 0x1f) {
  725.     int k1 = 32 - k;
  726.     unsigned long z = *--x;
  727.     if ((*x1 = (z >> k1)) != 0) {
  728.       ++n1;
  729.     }
  730.     while (x > xe) {
  731.       unsigned long w = *--x;
  732.       *--x1 = (z << k) | (w >> k1);
  733.       z = w;
  734.     }
  735.     *--x1 = z << k;
  736.   }
  737.   else
  738.     do {
  739.       *--x1 = *--x;
  740.     } while(x > xe);
  741.   while (x1 > xe)
  742.     *--x1 = 0;
  743.   b->wds = n1 - 1;
  744.   return b;
  745. }
  746.  
  747. static int
  748. cmp
  749. #ifdef KR_headers
  750.         (a, b) Bigint *a, *b;
  751. #else
  752.         (Bigint *a, Bigint *b)
  753. #endif
  754. {
  755.         unsigned long *xa, *xa0, *xb, *xb0;
  756.         int i, j;
  757.  
  758.         i = a->wds;
  759.         j = b->wds;
  760. #ifdef DEBUG
  761.         if (i > 1 && !a->x[i-1])
  762.                 Bug("cmp called with a->x[a->wds-1] == 0");
  763.         if (j > 1 && !b->x[j-1])
  764.                 Bug("cmp called with b->x[b->wds-1] == 0");
  765. #endif
  766.         if (i -= j)
  767.                 return i;
  768.         xa0 = a->x;
  769.         xa = xa0 + j;
  770.         xb0 = b->x;
  771.         xb = xb0 + j;
  772.         for(;;) {
  773.                 if (*--xa != *--xb)
  774.                         return *xa < *xb ? -1 : 1;
  775.                 if (xa <= xa0)
  776.                         break;
  777.                 }
  778.         return 0;
  779.         }
  780.  
  781. /* Do: c = a-b. */
  782.  
  783. static Bigint *
  784. diff
  785. #ifdef KR_headers
  786.         (c, a, b) Bigint *c, *a, *b;
  787. #else
  788.         (Bigint *c, Bigint *a, Bigint *b)
  789. #endif
  790. {
  791.         int i, wa, wb;
  792.         long borrow, y; /* We need signed shifts here. */
  793.         unsigned long *xa, *xae, *xb, *xbe, *xc;
  794.         long z;
  795.  
  796.         i = cmp(a,b);
  797.         if (!i) {
  798.                 c = Brealloc(c, 0);
  799.                 c->wds = 1;
  800.                 c->x[0] = 0;
  801.                 return c;
  802.                 }
  803.         if (i < 0) {
  804.                 Bigint *tmp = a;
  805.                 a = b;
  806.                 b = tmp;
  807.                 i = 1;
  808.                 }
  809.         else
  810.                 i = 0;
  811.         c = Brealloc(c, a->k);
  812.         c->sign = i;
  813.         wa = a->wds;
  814.         xa = a->x;
  815.         xae = xa + wa;
  816.         wb = b->wds;
  817.         xb = b->x;
  818.         xbe = xb + wb;
  819.         xc = c->x;
  820.         borrow = 0;
  821.         do {
  822.                 y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
  823.                 borrow = y >> 16;
  824.                 Sign_Extend(borrow, y);
  825.                 z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
  826.                 borrow = z >> 16;
  827.                 Sign_Extend(borrow, z);
  828.                 Storeinc(xc, z, y);
  829.                 }
  830.                 while(xb < xbe);
  831.         while(xa < xae) {
  832.                 y = (*xa & 0xffff) + borrow;
  833.                 borrow = y >> 16;
  834.                 Sign_Extend(borrow, y);
  835.                 z = (*xa++ >> 16) + borrow;
  836.                 borrow = z >> 16;
  837.                 Sign_Extend(borrow, z);
  838.                 Storeinc(xc, z, y);
  839.                 }
  840.         while(!*--xc)
  841.                 wa--;
  842.         c->wds = wa;
  843.         return c;
  844.         }
  845.  
  846. static double
  847. ulp
  848. #ifdef KR_headers
  849.         (x) double x;
  850. #else
  851.         (double x)
  852. #endif
  853. {
  854.         register long L;
  855.         double a;
  856.  
  857.         L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
  858. #ifndef Sudden_Underflow
  859.         if (L > 0) {
  860. #endif
  861. #ifdef IBM
  862.                 L |= Exp_msk1 >> 4;
  863. #endif
  864.                 word0(a) = L;
  865.                 word1(a) = 0;
  866. #ifndef Sudden_Underflow
  867.                 }
  868.         else {
  869.                 L = -L >> Exp_shift;
  870.                 if (L < Exp_shift) {
  871.                         word0(a) = 0x80000 >> L;
  872.                         word1(a) = 0;
  873.                         }
  874.                 else {
  875.                         word0(a) = 0;
  876.                         L -= Exp_shift;
  877.                         word1(a) = L >= 31 ? 1 : 1 << 31 - L;
  878.                         }
  879.                 }
  880. #endif
  881.         return a;
  882.         }
  883.  
  884. static double
  885. b2d
  886. #ifdef KR_headers
  887.         (a, e) Bigint *a; int *e;
  888. #else
  889.         (Bigint *a, int *e)
  890. #endif
  891. {
  892.         unsigned long *xa, *xa0, w, y, z;
  893.         int k;
  894.         double d;
  895. #ifdef VAX
  896.         unsigned long d0, d1;
  897. #else
  898. #define d0 word0(d)
  899. #define d1 word1(d)
  900. #endif
  901.  
  902.         xa0 = a->x;
  903.         xa = xa0 + a->wds;
  904.         y = *--xa;
  905. #ifdef DEBUG
  906.         if (!y) Bug("zero y in b2d");
  907. #endif
  908.         k = hi0bits(y);
  909.         *e = 32 - k;
  910.         if (k < Ebits) {
  911.                 d0 = Exp_1 | y >> Ebits - k;
  912.                 w = xa > xa0 ? *--xa : 0;
  913.                 d1 = y << (32-Ebits) + k | w >> Ebits - k;
  914.                 goto ret_d;
  915.                 }
  916.         z = xa > xa0 ? *--xa : 0;
  917.         if (k -= Ebits) {
  918.                 d0 = Exp_1 | y << k | z >> 32 - k;
  919.                 y = xa > xa0 ? *--xa : 0;
  920.                 d1 = z << k | y >> 32 - k;
  921.                 }
  922.         else {
  923.                 d0 = Exp_1 | y;
  924.                 d1 = z;
  925.                 }
  926.  ret_d:
  927. #ifdef VAX
  928.         word0(d) = d0 >> 16 | d0 << 16;
  929.         word1(d) = d1 >> 16 | d1 << 16;
  930. #else
  931. #undef d0
  932. #undef d1
  933. #endif
  934.         return d;
  935.         }
  936.  
  937. static Bigint *
  938. d2b
  939. #ifdef KR_headers
  940.         (result, d, e, bits) Bigint *result; double d; int *e, *bits;
  941. #else
  942.         (Bigint *result, double d, int *e, int *bits)
  943. #endif
  944. {
  945.         int de, i, k;
  946.         unsigned long *x, y, z;
  947. #ifdef VAX
  948.         unsigned long d0, d1;
  949.         d0 = word0(d) >> 16 | word0(d) << 16;
  950.         d1 = word1(d) >> 16 | word1(d) << 16;
  951. #else
  952. #define d0 word0(d)
  953. #define d1 word1(d)
  954. #endif
  955.  
  956.         result = Brealloc(result, 1);
  957.         x = result->x;
  958.  
  959.         z = d0 & Frac_mask;
  960.         d0 &= 0x7fffffff;       /* clear sign bit, which we ignore */
  961.  
  962.         de = (int)(d0 >> Exp_shift);  /* The exponent part of d. */
  963.  
  964.     /* Put back the suppressed high-order bit, if normalized. */
  965. #ifndef IBM
  966. #ifndef Sudden_Underflow
  967.         if (de)
  968. #endif
  969.       z |= Exp_msk11;
  970. #endif
  971.  
  972.         if (y = d1) {
  973.                 if (k = lo0bits(&y)) {
  974.                         x[0] = y | z << 32 - k;
  975.                         z >>= k;
  976.                         }
  977.                 else
  978.                         x[0] = y;
  979.                 i = result->wds = (x[1] = z) ? 2 : 1;
  980.                 }
  981.         else {
  982. #ifdef DEBUG
  983.                 if (!z)
  984.                         Bug("Zero passed to d2b");
  985. #endif
  986.                 k = lo0bits(&z);
  987.                 x[0] = z;
  988.                 i = result->wds = 1;
  989.                 k += 32;
  990.                 }
  991. #ifndef Sudden_Underflow
  992.         if (de) {
  993. #endif
  994. #ifdef IBM
  995.                 *e = (de - Bias - (P-1) << 2) + k;
  996.                 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
  997. #else
  998.                 *e = de - Bias - (P-1) + k;
  999.                 *bits = P - k;
  1000. #endif
  1001. #ifndef Sudden_Underflow
  1002.                 }
  1003.         else {
  1004.                 *e = de - Bias - (P-1) + 1 + k;
  1005.                 *bits = 32*i - hi0bits(x[i-1]);
  1006.                 }
  1007. #endif
  1008.         return result;
  1009.         }
  1010. #undef d0
  1011. #undef d1
  1012.  
  1013. static double
  1014. ratio
  1015. #ifdef KR_headers
  1016.         (a, b) Bigint *a, *b;
  1017. #else
  1018.         (Bigint *a, Bigint *b)
  1019. #endif
  1020. {
  1021.         double da, db;
  1022.         int k, ka, kb;
  1023.  
  1024.         da = b2d(a, &ka);
  1025.         db = b2d(b, &kb);
  1026.         k = ka - kb + 32*(a->wds - b->wds);
  1027. #ifdef IBM
  1028.         if (k > 0) {
  1029.                 word0(da) += (k >> 2)*Exp_msk1;
  1030.                 if (k &= 3)
  1031.                         da *= 1 << k;
  1032.                 }
  1033.         else {
  1034.                 k = -k;
  1035.                 word0(db) += (k >> 2)*Exp_msk1;
  1036.                 if (k &= 3)
  1037.                         db *= 1 << k;
  1038.                 }
  1039. #else
  1040.         if (k > 0)
  1041.                 word0(da) += k*Exp_msk1;
  1042.         else {
  1043.                 k = -k;
  1044.                 word0(db) += k*Exp_msk1;
  1045.                 }
  1046. #endif
  1047.         return da / db;
  1048.         }
  1049.  
  1050. static double
  1051. tens[] = {
  1052.                 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
  1053.                 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
  1054.                 1e20, 1e21, 1e22
  1055. #ifdef VAX
  1056.                 , 1e23, 1e24
  1057. #endif
  1058.                 };
  1059.  
  1060. static double
  1061. #ifdef IEEE_Arith
  1062. bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
  1063. static double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
  1064. #define n_bigtens 5
  1065. #else
  1066. #ifdef IBM
  1067. bigtens[] = { 1e16, 1e32, 1e64 };
  1068. static double tinytens[] = { 1e-16, 1e-32, 1e-64 };
  1069. #define n_bigtens 3
  1070. #else
  1071. bigtens[] = { 1e16, 1e32 };
  1072. static double tinytens[] = { 1e-16, 1e-32 };
  1073. #define n_bigtens 2
  1074. #endif
  1075. #endif
  1076.  
  1077.  double
  1078. _IO_strtod
  1079. #ifdef KR_headers
  1080.         (s00, se) CONST char *s00; char **se;
  1081. #else
  1082.         (CONST char *s00, char **se)
  1083. #endif
  1084. {
  1085.         int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
  1086.                  e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
  1087.         CONST char *s, *s0, *s1;
  1088.         double aadj, aadj1, adj, rv, rv0;
  1089.         long L;
  1090.         unsigned long y, z;
  1091.     Bigint _bb, _b_avail, _bd, _bd0, _bs, _delta;
  1092.     Bigint *bb = Binit(&_bb);
  1093.     Bigint *bd = Binit(&_bd);
  1094.     Bigint *bd0 = Binit(&_bd0);
  1095.     Bigint *bs = Binit(&_bs);
  1096.     Bigint *b_avail = Binit(&_b_avail);
  1097.     Bigint *delta = Binit(&_delta);
  1098.  
  1099.     TEST_ENDIANNESS;
  1100.         sign = nz0 = nz = 0;
  1101.         rv = 0.;
  1102.         for(s = s00;;s++) switch(*s) {
  1103.                 case '-':
  1104.                         sign = 1;
  1105.                         /* no break */
  1106.                 case '+':
  1107.                         if (*++s)
  1108.                                 goto break2;
  1109.                         /* no break */
  1110.                 case 0:
  1111.                         goto ret;
  1112.                 case '\t':
  1113.                 case '\n':
  1114.                 case '\v':
  1115.                 case '\f':
  1116.                 case '\r':
  1117.                 case ' ':
  1118.                         continue;
  1119.                 default:
  1120.                         goto break2;
  1121.                 }
  1122.  break2:
  1123.         if (*s == '0') {
  1124.                 nz0 = 1;
  1125.                 while(*++s == '0') ;
  1126.                 if (!*s)
  1127.                         goto ret;
  1128.                 }
  1129.         s0 = s;
  1130.         y = z = 0;
  1131.         for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
  1132.                 if (nd < 9)
  1133.                         y = 10*y + c - '0';
  1134.                 else if (nd < 16)
  1135.                         z = 10*z + c - '0';
  1136.         nd0 = nd;
  1137.         if (c == '.') {
  1138.                 c = *++s;
  1139.                 if (!nd) {
  1140.                         for(; c == '0'; c = *++s)
  1141.                                 nz++;
  1142.                         if (c > '0' && c <= '9') {
  1143.                                 s0 = s;
  1144.                                 nf += nz;
  1145.                                 nz = 0;
  1146.                                 goto have_dig;
  1147.                                 }
  1148.                         goto dig_done;
  1149.                         }
  1150.                 for(; c >= '0' && c <= '9'; c = *++s) {
  1151.  have_dig:
  1152.                         nz++;
  1153.                         if (c -= '0') {
  1154.                                 nf += nz;
  1155.                                 for(i = 1; i < nz; i++)
  1156.                                         if (nd++ < 9)
  1157.                                                 y *= 10;
  1158.                                         else if (nd <= DBL_DIG + 1)
  1159.                                                 z *= 10;
  1160.                                 if (nd++ < 9)
  1161.                                         y = 10*y + c;
  1162.                                 else if (nd <= DBL_DIG + 1)
  1163.                                         z = 10*z + c;
  1164.                                 nz = 0;
  1165.                                 }
  1166.                         }
  1167.                 }
  1168.  dig_done:
  1169.         e = 0;
  1170.         if (c == 'e' || c == 'E') {
  1171.                 if (!nd && !nz && !nz0) {
  1172.                         s = s00;
  1173.                         goto ret;
  1174.                         }
  1175.                 s00 = s;
  1176.                 esign = 0;
  1177.                 switch(c = *++s) {
  1178.                         case '-':
  1179.                                 esign = 1;
  1180.                         case '+':
  1181.                                 c = *++s;
  1182.                         }
  1183.                 if (c >= '0' && c <= '9') {
  1184.                         while(c == '0')
  1185.                                 c = *++s;
  1186.                         if (c > '0' && c <= '9') {
  1187.                                 e = c - '0';
  1188.                                 s1 = s;
  1189.                                 while((c = *++s) >= '0' && c <= '9')
  1190.                                         e = 10*e + c - '0';
  1191.                                 if (s - s1 > 8)
  1192.                                         /* Avoid confusion from exponents
  1193.                                          * so large that e might overflow.
  1194.                                          */
  1195.                                         e = 9999999;
  1196.                                 if (esign)
  1197.                                         e = -e;
  1198.                                 }
  1199.                         else
  1200.                                 e = 0;
  1201.                         }
  1202.                 else
  1203.                         s = s00;
  1204.                 }
  1205.         if (!nd) {
  1206.                 if (!nz && !nz0)
  1207.                         s = s00;
  1208.                 goto ret;
  1209.                 }
  1210.         e1 = e -= nf;
  1211.  
  1212.         /* Now we have nd0 digits, starting at s0, followed by a
  1213.          * decimal point, followed by nd-nd0 digits.  The number we're
  1214.          * after is the integer represented by those digits times
  1215.          * 10**e */
  1216.  
  1217.         if (!nd0)
  1218.                 nd0 = nd;
  1219.         k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
  1220.         rv = y;
  1221.         if (k > 9)
  1222.                 rv = tens[k - 9] * rv + z;
  1223.         if (nd <= DBL_DIG
  1224. #ifndef RND_PRODQUOT
  1225.                 && FLT_ROUNDS == 1
  1226. #endif
  1227.                         ) {
  1228.                 if (!e)
  1229.                         goto ret;
  1230.                 if (e > 0) {
  1231.                         if (e <= Ten_pmax) {
  1232. #ifdef VAX
  1233.                                 goto vax_ovfl_check;
  1234. #else
  1235.                                 /* rv = */ rounded_product(rv, tens[e]);
  1236.                                 goto ret;
  1237. #endif
  1238.                                 }
  1239.                         i = DBL_DIG - nd;
  1240.                         if (e <= Ten_pmax + i) {
  1241.                                 /* A fancier test would sometimes let us do
  1242.                                  * this for larger i values.
  1243.                                  */
  1244.                                 e -= i;
  1245.                                 rv *= tens[i];
  1246. #ifdef VAX
  1247.                                 /* VAX exponent range is so narrow we must
  1248.                                  * worry about overflow here...
  1249.                                  */
  1250.  vax_ovfl_check:
  1251.                                 word0(rv) -= P*Exp_msk1;
  1252.                                 /* rv = */ rounded_product(rv, tens[e]);
  1253.                                 if ((word0(rv) & Exp_mask)
  1254.                                  > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
  1255.                                         goto ovfl;
  1256.                                 word0(rv) += P*Exp_msk1;
  1257. #else
  1258.                                 /* rv = */ rounded_product(rv, tens[e]);
  1259. #endif
  1260.                                 goto ret;
  1261.                                 }
  1262.                         }
  1263. #ifndef Inaccurate_Divide
  1264.                 else if (e >= -Ten_pmax) {
  1265.                         /* rv = */ rounded_quotient(rv, tens[-e]);
  1266.                         goto ret;
  1267.                         }
  1268. #endif
  1269.                 }
  1270.         e1 += nd - k;
  1271.  
  1272.         /* Get starting approximation = rv * 10**e1 */
  1273.  
  1274.         if (e1 > 0) {
  1275.                 if (i = e1 & 15)
  1276.                         rv *= tens[i];
  1277.                 if (e1 &= ~15) {
  1278.                         if (e1 > DBL_MAX_10_EXP) {
  1279.  ovfl:
  1280.                                 errno = ERANGE;
  1281. #if defined(sun) && !defined(__svr4__)
  1282. /* SunOS defines HUGE_VAL as __infinity(), which is in libm. */
  1283. #undef HUGE_VAL
  1284. #endif
  1285. #ifndef HUGE_VAL
  1286. #define HUGE_VAL        1.7976931348623157E+308
  1287. #endif
  1288.                                 rv = HUGE_VAL;
  1289.                                 goto ret;
  1290.                                 }
  1291.                         if (e1 >>= 4) {
  1292.                                 for(j = 0; e1 > 1; j++, e1 >>= 1)
  1293.                                         if (e1 & 1)
  1294.                                                 rv *= bigtens[j];
  1295.                         /* The last multiplication could overflow. */
  1296.                                 word0(rv) -= P*Exp_msk1;
  1297.                                 rv *= bigtens[j];
  1298.                                 if ((z = word0(rv) & Exp_mask)
  1299.                                  > Exp_msk1*(DBL_MAX_EXP+Bias-P))
  1300.                                         goto ovfl;
  1301.                                 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
  1302.                                         /* set to largest number */
  1303.                                         /* (Can't trust DBL_MAX) */
  1304.                                         word0(rv) = Big0;
  1305.                                         word1(rv) = Big1;
  1306.                                         }
  1307.                                 else
  1308.                                         word0(rv) += P*Exp_msk1;
  1309.                                 }
  1310.  
  1311.                         }
  1312.                 }
  1313.         else if (e1 < 0) {
  1314.                 e1 = -e1;
  1315.                 if (i = e1 & 15)
  1316.                         rv /= tens[i];
  1317.                 if (e1 &= ~15) {
  1318.                         e1 >>= 4;
  1319.                         for(j = 0; e1 > 1; j++, e1 >>= 1)
  1320.                                 if (e1 & 1)
  1321.                                         rv *= tinytens[j];
  1322.                         /* The last multiplication could underflow. */
  1323.                         rv0 = rv;
  1324.                         rv *= tinytens[j];
  1325.                         if (!rv) {
  1326.                                 rv = 2.*rv0;
  1327.                                 rv *= tinytens[j];
  1328.                                 if (!rv) {
  1329.  undfl:
  1330.                                         rv = 0.;
  1331.                                         errno = ERANGE;
  1332.                                         goto ret;
  1333.                                         }
  1334.                                 word0(rv) = Tiny0;
  1335.                                 word1(rv) = Tiny1;
  1336.                                 /* The refinement below will clean
  1337.                                  * this approximation up.
  1338.                                  */
  1339.                                 }
  1340.                         }
  1341.                 }
  1342.  
  1343.         /* Now the hard part -- adjusting rv to the correct value.*/
  1344.  
  1345.         /* Put digits into bd: true value = bd * 10^e */
  1346.  
  1347.         bd0 = s2b(bd0, s0, nd0, nd, y);
  1348.     bd = Brealloc(bd, bd0->k);
  1349.  
  1350.         for(;;) {
  1351.                 Bcopy(bd, bd0);
  1352.                 bb = d2b(bb, rv, &bbe, &bbbits);    /* rv = bb * 2^bbe */
  1353.                 bs = i2b(bs, 1);
  1354.  
  1355.                 if (e >= 0) {
  1356.                         bb2 = bb5 = 0;
  1357.                         bd2 = bd5 = e;
  1358.                         }
  1359.                 else {
  1360.                         bb2 = bb5 = -e;
  1361.                         bd2 = bd5 = 0;
  1362.                         }
  1363.                 if (bbe >= 0)
  1364.                         bb2 += bbe;
  1365.                 else
  1366.                         bd2 -= bbe;
  1367.                 bs2 = bb2;
  1368. #ifdef Sudden_Underflow
  1369. #ifdef IBM
  1370.                 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
  1371. #else
  1372.                 j = P + 1 - bbbits;
  1373. #endif
  1374. #else
  1375.                 i = bbe + bbbits - 1;   /* logb(rv) */
  1376.                 if (i < Emin)   /* denormal */
  1377.                         j = bbe + (P-Emin);
  1378.                 else
  1379.                         j = P + 1 - bbbits;
  1380. #endif
  1381.                 bb2 += j;
  1382.                 bd2 += j;
  1383.                 i = bb2 < bd2 ? bb2 : bd2;
  1384.                 if (i > bs2)
  1385.                         i = bs2;
  1386.                 if (i > 0) {
  1387.                         bb2 -= i;
  1388.                         bd2 -= i;
  1389.                         bs2 -= i;
  1390.                         }
  1391.                 if (bb5 > 0) {
  1392.             Bigint *b_tmp;
  1393.                         bs = pow5mult(bs, bb5);
  1394.                         b_tmp = mult(b_avail, bs, bb);
  1395.                         b_avail = bb;
  1396.                         bb = b_tmp;
  1397.                         }
  1398.                 if (bb2 > 0)
  1399.                         bb = lshift(bb, bb2);
  1400.                 if (bd5 > 0)
  1401.                         bd = pow5mult(bd, bd5);
  1402.                 if (bd2 > 0)
  1403.                         bd = lshift(bd, bd2);
  1404.                 if (bs2 > 0)
  1405.                         bs = lshift(bs, bs2);
  1406.                 delta = diff(delta, bb, bd);
  1407.                 dsign = delta->sign;
  1408.                 delta->sign = 0;
  1409.                 i = cmp(delta, bs);
  1410.                 if (i < 0) {
  1411.                         /* Error is less than half an ulp -- check for
  1412.                          * special case of mantissa a power of two.
  1413.                          */
  1414.                         if (dsign || word1(rv) || word0(rv) & Bndry_mask)
  1415.                                 break;
  1416.                         delta = lshift(delta,Log2P);
  1417.                         if (cmp(delta, bs) > 0)
  1418.                                 goto drop_down;
  1419.                         break;
  1420.                         }
  1421.                 if (i == 0) {
  1422.                         /* exactly half-way between */
  1423.                         if (dsign) {
  1424.                                 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
  1425.                                  &&  word1(rv) == 0xffffffff) {
  1426.                                         /*boundary case -- increment exponent*/
  1427.                                         word0(rv) = (word0(rv) & Exp_mask)
  1428.                                                 + Exp_msk1
  1429. #ifdef IBM
  1430.                                                 | Exp_msk1 >> 4
  1431. #endif
  1432.                                                 ;
  1433.                                         word1(rv) = 0;
  1434.                                         break;
  1435.                                         }
  1436.                                 }
  1437.                         else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
  1438.  drop_down:
  1439.                                 /* boundary case -- decrement exponent */
  1440. #ifdef Sudden_Underflow
  1441.                                 L = word0(rv) & Exp_mask;
  1442. #ifdef IBM
  1443.                                 if (L <  Exp_msk1)
  1444. #else
  1445.                                 if (L <= Exp_msk1)
  1446. #endif
  1447.                                         goto undfl;
  1448.                                 L -= Exp_msk1;
  1449. #else
  1450.                                 L = (word0(rv) & Exp_mask) - Exp_msk1;
  1451. #endif
  1452.                                 word0(rv) = L | Bndry_mask1;
  1453.                                 word1(rv) = 0xffffffff;
  1454. #ifdef IBM
  1455.                                 continue;
  1456. #else
  1457.                                 break;
  1458. #endif
  1459.                                 }
  1460. #ifndef ROUND_BIASED
  1461.                         if (!(word1(rv) & LSB))
  1462.                                 break;
  1463. #endif
  1464.                         if (dsign)
  1465.                                 rv += ulp(rv);
  1466. #ifndef ROUND_BIASED
  1467.                         else {
  1468.                                 rv -= ulp(rv);
  1469. #ifndef Sudden_Underflow
  1470.                                 if (!rv)
  1471.                                         goto undfl;
  1472. #endif
  1473.                                 }
  1474. #endif
  1475.                         break;
  1476.                         }
  1477.                 if ((aadj = ratio(delta, bs)) <= 2.) {
  1478.                         if (dsign)
  1479.                                 aadj = aadj1 = 1.;
  1480.                         else if (word1(rv) || word0(rv) & Bndry_mask) {
  1481. #ifndef Sudden_Underflow
  1482.                                 if (word1(rv) == Tiny1 && !word0(rv))
  1483.                                         goto undfl;
  1484. #endif
  1485.                                 aadj = 1.;
  1486.                                 aadj1 = -1.;
  1487.                                 }
  1488.                         else {
  1489.                                 /* special case -- power of FLT_RADIX to be */
  1490.                                 /* rounded down... */
  1491.  
  1492.                                 if (aadj < 2./FLT_RADIX)
  1493.                                         aadj = 1./FLT_RADIX;
  1494.                                 else
  1495.                                         aadj *= 0.5;
  1496.                                 aadj1 = -aadj;
  1497.                                 }
  1498.                         }
  1499.                 else {
  1500.                         aadj *= 0.5;
  1501.                         aadj1 = dsign ? aadj : -aadj;
  1502. #ifdef Check_FLT_ROUNDS
  1503.                         switch(FLT_ROUNDS) {
  1504.                                 case 2: /* towards +infinity */
  1505.                                         aadj1 -= 0.5;
  1506.                                         break;
  1507.                                 case 0: /* towards 0 */
  1508.                                 case 3: /* towards -infinity */
  1509.                                         aadj1 += 0.5;
  1510.                                 }
  1511. #else
  1512.                         if (FLT_ROUNDS == 0)
  1513.                                 aadj1 += 0.5;
  1514. #endif
  1515.                         }
  1516.                 y = word0(rv) & Exp_mask;
  1517.  
  1518.                 /* Check for overflow */
  1519.  
  1520.                 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
  1521.                         rv0 = rv;
  1522.                         word0(rv) -= P*Exp_msk1;
  1523.                         adj = aadj1 * ulp(rv);
  1524.                         rv += adj;
  1525.                         if ((word0(rv) & Exp_mask) >=
  1526.                                         Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
  1527.                                 if (word0(rv0) == Big0 && word1(rv0) == Big1)
  1528.                                         goto ovfl;
  1529.                                 word0(rv) = Big0;
  1530.                                 word1(rv) = Big1;
  1531.                                 continue;
  1532.                                 }
  1533.                         else
  1534.                                 word0(rv) += P*Exp_msk1;
  1535.                         }
  1536.                 else {
  1537. #ifdef Sudden_Underflow
  1538.                         if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
  1539.                                 rv0 = rv;
  1540.                                 word0(rv) += P*Exp_msk1;
  1541.                                 adj = aadj1 * ulp(rv);
  1542.                                 rv += adj;
  1543. #ifdef IBM
  1544.                                 if ((word0(rv) & Exp_mask) <  P*Exp_msk1)
  1545. #else
  1546.                                 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
  1547. #endif
  1548.                                         {
  1549.                                         if (word0(rv0) == Tiny0
  1550.                                          && word1(rv0) == Tiny1)
  1551.                                                 goto undfl;
  1552.                                         word0(rv) = Tiny0;
  1553.                                         word1(rv) = Tiny1;
  1554.                                         continue;
  1555.                                         }
  1556.                                 else
  1557.                                         word0(rv) -= P*Exp_msk1;
  1558.                                 }
  1559.                         else {
  1560.                                 adj = aadj1 * ulp(rv);
  1561.                                 rv += adj;
  1562.                                 }
  1563. #else
  1564.                         /* Compute adj so that the IEEE rounding rules will
  1565.                          * correctly round rv + adj in some half-way cases.
  1566.                          * If rv * ulp(rv) is denormalized (i.e.,
  1567.                          * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
  1568.                          * trouble from bits lost to denormalization;
  1569.                          * example: 1.2e-307 .
  1570.                          */
  1571.                         if (y <= (P-1)*Exp_msk1 && aadj >= 1.) {
  1572.                                 aadj1 = (double)(int)(aadj + 0.5);
  1573.                                 if (!dsign)
  1574.                                         aadj1 = -aadj1;
  1575.                                 }
  1576.                         adj = aadj1 * ulp(rv);
  1577.                         rv += adj;
  1578. #endif
  1579.                         }
  1580.                 z = word0(rv) & Exp_mask;
  1581.                 if (y == z) {
  1582.                         /* Can we stop now? */
  1583.                         L = (long)aadj;
  1584.                         aadj -= L;
  1585.                         /* The tolerances below are conservative. */
  1586.                         if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
  1587.                                 if (aadj < .4999999 || aadj > .5000001)
  1588.                                         break;
  1589.                                 }
  1590.                         else if (aadj < .4999999/FLT_RADIX)
  1591.                                 break;
  1592.                         }
  1593.                 }
  1594.         Bfree(bb);
  1595.         Bfree(bd);
  1596.         Bfree(bs);
  1597.         Bfree(bd0);
  1598.         Bfree(delta);
  1599.     Bfree(b_avail);
  1600.  ret:
  1601.         if (se)
  1602.                 *se = (char *)s;
  1603.         return sign ? -rv : rv;
  1604.         }
  1605.  
  1606. static int
  1607. quorem
  1608. #ifdef KR_headers
  1609.         (b, S) Bigint *b, *S;
  1610. #else
  1611.         (Bigint *b, Bigint *S)
  1612. #endif
  1613. {
  1614.         int n;
  1615.         long borrow, y;
  1616.         unsigned long carry, q, ys;
  1617.         unsigned long *bx, *bxe, *sx, *sxe;
  1618.         long z;
  1619.         unsigned long si, zs;
  1620.  
  1621.         n = S->wds;
  1622. #ifdef DEBUG
  1623.         /*debug*/ if (b->wds > n)
  1624.         /*debug*/       Bug("oversize b in quorem");
  1625. #endif
  1626.         if (b->wds < n)
  1627.                 return 0;
  1628.         sx = S->x;
  1629.         sxe = sx + --n;
  1630.         bx = b->x;
  1631.         bxe = bx + n;
  1632.         q = *bxe / (*sxe + 1);  /* ensure q <= true quotient */
  1633. #ifdef DEBUG
  1634.         /*debug*/ if (q > 9)
  1635.         /*debug*/       Bug("oversized quotient in quorem");
  1636. #endif
  1637.         if (q) {
  1638.                 borrow = 0;
  1639.                 carry = 0;
  1640.                 do {
  1641.                         si = *sx++;
  1642.                         ys = (si & 0xffff) * q + carry;
  1643.                         zs = (si >> 16) * q + (ys >> 16);
  1644.                         carry = zs >> 16;
  1645.                         y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
  1646.                         borrow = y >> 16;
  1647.                         Sign_Extend(borrow, y);
  1648.                         z = (*bx >> 16) - (zs & 0xffff) + borrow;
  1649.                         borrow = z >> 16;
  1650.                         Sign_Extend(borrow, z);
  1651.                         Storeinc(bx, z, y);
  1652.                         }
  1653.                         while(sx <= sxe);
  1654.                 if (!*bxe) {
  1655.                         bx = b->x;
  1656.                         while(--bxe > bx && !*bxe)
  1657.                                 --n;
  1658.                         b->wds = n;
  1659.                         }
  1660.                 }
  1661.         if (cmp(b, S) >= 0) {
  1662.                 q++;
  1663.                 borrow = 0;
  1664.                 carry = 0;
  1665.                 bx = b->x;
  1666.                 sx = S->x;
  1667.                 do {
  1668.                         si = *sx++;
  1669.                         ys = (si & 0xffff) + carry;
  1670.                         zs = (si >> 16) + (ys >> 16);
  1671.                         carry = zs >> 16;
  1672.                         y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
  1673.                         borrow = y >> 16;
  1674.                         Sign_Extend(borrow, y);
  1675.                         z = (*bx >> 16) - (zs & 0xffff) + borrow;
  1676.                         borrow = z >> 16;
  1677.                         Sign_Extend(borrow, z);
  1678.                         Storeinc(bx, z, y);
  1679.                         }
  1680.                         while(sx <= sxe);
  1681.                 bx = b->x;
  1682.                 bxe = bx + n;
  1683.                 if (!*bxe) {
  1684.                         while(--bxe > bx && !*bxe)
  1685.                                 --n;
  1686.                         b->wds = n;
  1687.                         }
  1688.                 }
  1689.         return q;
  1690.         }
  1691.  
  1692. /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
  1693.  *
  1694.  * Inspired by "How to Print Floating-Point Numbers Accurately" by
  1695.  * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
  1696.  *
  1697.  * Modifications:
  1698.  *      1. Rather than iterating, we use a simple numeric overestimate
  1699.  *         to determine k = floor(log10(d)).  We scale relevant
  1700.  *         quantities using O(log2(k)) rather than O(k) multiplications.
  1701.  *      2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
  1702.  *         try to generate digits strictly left to right.  Instead, we
  1703.  *         compute with fewer bits and propagate the carry if necessary
  1704.  *         when rounding the final digit up.  This is often faster.
  1705.  *      3. Under the assumption that input will be rounded nearest,
  1706.  *         mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
  1707.  *         That is, we allow equality in stopping tests when the
  1708.  *         round-nearest rule will give the same floating-point value
  1709.  *         as would satisfaction of the stopping test with strict
  1710.  *         inequality.
  1711.  *      4. We remove common factors of powers of 2 from relevant
  1712.  *         quantities.
  1713.  *      5. When converting floating-point integers less than 1e16,
  1714.  *         we use floating-point arithmetic rather than resorting
  1715.  *         to multiple-precision integers.
  1716.  *      6. When asked to produce fewer than 15 digits, we first try
  1717.  *         to get by with floating-point arithmetic; we resort to
  1718.  *         multiple-precision integer arithmetic only if we cannot
  1719.  *         guarantee that the floating-point calculation has given
  1720.  *         the correctly rounded result.  For k requested digits and
  1721.  *         "uniformly" distributed input, the probability is
  1722.  *         something like 10^(k-15) that we must resort to the long
  1723.  *         calculation.
  1724.  */
  1725.  
  1726.  char *
  1727. _IO_dtoa
  1728. #ifdef KR_headers
  1729.         (d, mode, ndigits, decpt, sign, rve)
  1730.         double d; int mode, ndigits, *decpt, *sign; char **rve;
  1731. #else
  1732.         (double d, int mode, int ndigits, int *decpt, int *sign, char **rve)
  1733. #endif
  1734. {
  1735.  /*     Arguments ndigits, decpt, sign are similar to those
  1736.         of ecvt and fcvt; trailing zeros are suppressed from
  1737.         the returned string.  If not null, *rve is set to point
  1738.         to the end of the return value.  If d is +-Infinity or NaN,
  1739.         then *decpt is set to 9999.
  1740.  
  1741.         mode:
  1742.                 0 ==> shortest string that yields d when read in
  1743.                         and rounded to nearest.
  1744.                 1 ==> like 0, but with Steele & White stopping rule;
  1745.                         e.g. with IEEE P754 arithmetic , mode 0 gives
  1746.                         1e23 whereas mode 1 gives 9.999999999999999e22.
  1747.                 2 ==> max(1,ndigits) significant digits.  This gives a
  1748.                         return value similar to that of ecvt, except
  1749.                         that trailing zeros are suppressed.
  1750.                 3 ==> through ndigits past the decimal point.  This
  1751.                         gives a return value similar to that from fcvt,
  1752.                         except that trailing zeros are suppressed, and
  1753.                         ndigits can be negative.
  1754.                 4-9 should give the same return values as 2-3, i.e.,
  1755.                         4 <= mode <= 9 ==> same return as mode
  1756.                         2 + (mode & 1).  These modes are mainly for
  1757.                         debugging; often they run slower but sometimes
  1758.                         faster than modes 2-3.
  1759.                 4,5,8,9 ==> left-to-right digit generation.
  1760.                 6-9 ==> don't try fast floating-point estimate
  1761.                         (if applicable).
  1762.  
  1763.                 Values of mode other than 0-9 are treated as mode 0.
  1764.  
  1765.                 Sufficient space is allocated to the return value
  1766.                 to hold the suppressed trailing zeros.
  1767.         */
  1768.  
  1769.         int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
  1770.                 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
  1771.                 spec_case, try_quick;
  1772.         long L;
  1773. #ifndef Sudden_Underflow
  1774.         int denorm;
  1775. #endif
  1776.     Bigint _b_avail, _b, _mhi, _mlo, _S;
  1777.     Bigint *b_avail = Binit(&_b_avail);
  1778.     Bigint *b = Binit(&_b);
  1779.     Bigint *S = Binit(&_S);
  1780.     /* mhi and mlo are only set and used if leftright. */
  1781.         Bigint *mhi = NULL, *mlo = NULL;
  1782.         double d2, ds, eps;
  1783.         char *s, *s0;
  1784.         static Bigint *result = NULL;
  1785.         static int result_k;
  1786.  
  1787.     TEST_ENDIANNESS;
  1788.         if (result) {
  1789.                 result->k = result_k;
  1790.                 result->maxwds = 1 << result_k;
  1791.                 }
  1792.  
  1793.         if (word0(d) & Sign_bit) {
  1794.                 /* set sign for everything, including 0's and NaNs */
  1795.                 *sign = 1;
  1796.                 word0(d) &= ~Sign_bit;  /* clear sign bit */
  1797.                 }
  1798.         else
  1799.                 *sign = 0;
  1800.  
  1801. #if defined(IEEE_Arith) + defined(VAX)
  1802. #ifdef IEEE_Arith
  1803.         if ((word0(d) & Exp_mask) == Exp_mask)
  1804. #else
  1805.         if (word0(d)  == 0x8000)
  1806. #endif
  1807.                 {
  1808.                 /* Infinity or NaN */
  1809.                 *decpt = 9999;
  1810. #ifdef IEEE_Arith
  1811.         if (!word1(d) && !(word0(d) & 0xfffff))
  1812.           {
  1813.             s = "Infinity";
  1814.             if (rve)
  1815.               *rve = s + 8;
  1816.           }
  1817.         else
  1818. #endif
  1819.           {
  1820.             s = "NaN";
  1821.             if (rve)
  1822.               *rve = s +3;
  1823.           }
  1824.                 return s;
  1825.                 }
  1826. #endif
  1827. #ifdef IBM
  1828.         d += 0; /* normalize */
  1829. #endif
  1830.         if (!d) {
  1831.                 *decpt = 1;
  1832.                 s = "0";
  1833.                 if (rve)
  1834.                         *rve = s + 1;
  1835.                 return s;
  1836.                 }
  1837.  
  1838.         b = d2b(b, d, &be, &bbits);
  1839.         i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
  1840. #ifndef Sudden_Underflow
  1841.         if (i) {
  1842. #endif
  1843.                 d2 = d;
  1844.                 word0(d2) &= Frac_mask1;
  1845.                 word0(d2) |= Exp_11;
  1846. #ifdef IBM
  1847.                 if (j = 11 - hi0bits(word0(d2) & Frac_mask))
  1848.                         d2 /= 1 << j;
  1849. #endif
  1850.  
  1851.                 i -= Bias;
  1852. #ifdef IBM
  1853.                 i <<= 2;
  1854.                 i += j;
  1855. #endif
  1856. #ifndef Sudden_Underflow
  1857.                 denorm = 0;
  1858.                 }
  1859.         else {
  1860.                 /* d is denormalized */
  1861.         unsigned long x;
  1862.  
  1863.                 i = bbits + be + (Bias + (P-1) - 1);
  1864.                 x = i > 32  ? word0(d) << 64 - i | word1(d) >> i - 32
  1865.                             : word1(d) << 32 - i;
  1866.                 d2 = x;
  1867.                 word0(d2) -= 31*Exp_msk1; /* adjust exponent */
  1868.                 i -= (Bias + (P-1) - 1) + 1;
  1869.                 denorm = 1;
  1870.                 }
  1871. #endif
  1872.  
  1873.     /* Now i is the unbiased base-2 exponent. */
  1874.  
  1875.         /* log(x)       ~=~ log(1.5) + (x-1.5)/1.5
  1876.          * log10(x)      =  log(x) / log(10)
  1877.          *              ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
  1878.          * log10(d) = i*log(2)/log(10) + log10(d2)
  1879.          *
  1880.          * This suggests computing an approximation k to log10(d) by
  1881.          *
  1882.          * k = i*0.301029995663981
  1883.          *      + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
  1884.          *
  1885.          * We want k to be too large rather than too small.
  1886.          * The error in the first-order Taylor series approximation
  1887.          * is in our favor, so we just round up the constant enough
  1888.          * to compensate for any error in the multiplication of
  1889.          * (i) by 0.301029995663981; since |i| <= 1077,
  1890.          * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
  1891.          * adding 1e-13 to the constant term more than suffices.
  1892.          * Hence we adjust the constant term to 0.1760912590558.
  1893.          * (We could get a more accurate k by invoking log10,
  1894.          *  but this is probably not worthwhile.)
  1895.          */
  1896.  
  1897.         ds = (d2-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
  1898.         k = (int)ds;
  1899.         if (ds < 0. && ds != k)
  1900.                 k--;    /* want k = floor(ds) */
  1901.         k_check = 1;
  1902.         if (k >= 0 && k <= Ten_pmax) {
  1903.                 if (d < tens[k])
  1904.                         k--;
  1905.                 k_check = 0;
  1906.                 }
  1907.         j = bbits - i - 1;
  1908.         if (j >= 0) {
  1909.                 b2 = 0;
  1910.                 s2 = j;
  1911.                 }
  1912.         else {
  1913.                 b2 = -j;
  1914.                 s2 = 0;
  1915.                 }
  1916.         if (k >= 0) {
  1917.                 b5 = 0;
  1918.                 s5 = k;
  1919.                 s2 += k;
  1920.                 }
  1921.         else {
  1922.                 b2 -= k;
  1923.                 b5 = -k;
  1924.                 s5 = 0;
  1925.                 }
  1926.         if (mode < 0 || mode > 9)
  1927.                 mode = 0;
  1928.         try_quick = 1;
  1929.         if (mode > 5) {
  1930.                 mode -= 4;
  1931.                 try_quick = 0;
  1932.                 }
  1933.         leftright = 1;
  1934.         switch(mode) {
  1935.                 case 0:
  1936.                 case 1:
  1937.                         ilim = ilim1 = -1;
  1938.                         i = 18;
  1939.                         ndigits = 0;
  1940.                         break;
  1941.                 case 2:
  1942.                         leftright = 0;
  1943.                         /* no break */
  1944.                 case 4:
  1945.                         if (ndigits <= 0)
  1946.                                 ndigits = 1;
  1947.                         ilim = ilim1 = i = ndigits;
  1948.                         break;
  1949.                 case 3:
  1950.                         leftright = 0;
  1951.                         /* no break */
  1952.                 case 5:
  1953.                         i = ndigits + k + 1;
  1954.                         ilim = i;
  1955.                         ilim1 = i - 1;
  1956.                         if (i <= 0)
  1957.                                 i = 1;
  1958.                 }
  1959.     /* i is now an upper bound of the number of digits to generate. */
  1960.         j = sizeof(unsigned long) * (1<<BIGINT_MINIMUM_K);
  1961.     /* The test is <= so as to allow room for the final '\0'. */
  1962.         for(result_k = BIGINT_MINIMUM_K; BIGINT_HEADER_SIZE + j <= i;
  1963.                 j <<= 1) result_k++;
  1964.         result = Brealloc(result, result_k);
  1965.         s = s0 = (char *)result;
  1966.  
  1967.         if (ilim >= 0 && ilim <= Quick_max && try_quick) {
  1968.  
  1969.                 /* Try to get by with floating-point arithmetic. */
  1970.  
  1971.                 i = 0;
  1972.                 d2 = d;
  1973.                 k0 = k;
  1974.                 ilim0 = ilim;
  1975.                 ieps = 2; /* conservative */
  1976.                 if (k > 0) {
  1977.                         ds = tens[k&0xf];
  1978.                         j = k >> 4;
  1979.                         if (j & Bletch) {
  1980.                                 /* prevent overflows */
  1981.                                 j &= Bletch - 1;
  1982.                                 d /= bigtens[n_bigtens-1];
  1983.                                 ieps++;
  1984.                                 }
  1985.                         for(; j; j >>= 1, i++)
  1986.                                 if (j & 1) {
  1987.                                         ieps++;
  1988.                                         ds *= bigtens[i];
  1989.                                         }
  1990.                         d /= ds;
  1991.                         }
  1992.                 else if (j1 = -k) {
  1993.                         d *= tens[j1 & 0xf];
  1994.                         for(j = j1 >> 4; j; j >>= 1, i++)
  1995.                                 if (j & 1) {
  1996.                                         ieps++;
  1997.                                         d *= bigtens[i];
  1998.                                         }
  1999.                         }
  2000.                 if (k_check && d < 1. && ilim > 0) {
  2001.                         if (ilim1 <= 0)
  2002.                                 goto fast_failed;
  2003.                         ilim = ilim1;
  2004.                         k--;
  2005.                         d *= 10.;
  2006.                         ieps++;
  2007.                         }
  2008.                 eps = ieps*d + 7.;
  2009.                 word0(eps) -= (P-1)*Exp_msk1;
  2010.                 if (ilim == 0) {
  2011.                         d -= 5.;
  2012.                         if (d > eps)
  2013.                                 goto one_digit;
  2014.                         if (d < -eps)
  2015.                                 goto no_digits;
  2016.                         goto fast_failed;
  2017.                         }
  2018. #ifndef No_leftright
  2019.                 if (leftright) {
  2020.                         /* Use Steele & White method of only
  2021.                          * generating digits needed.
  2022.                          */
  2023.                         eps = 0.5/tens[ilim-1] - eps;
  2024.                         for(i = 0;;) {
  2025.                                 L = (long)d;
  2026.                                 d -= L;
  2027.                                 *s++ = '0' + (int)L;
  2028.                                 if (d < eps)
  2029.                                         goto ret1;
  2030.                                 if (1. - d < eps)
  2031.                                         goto bump_up;
  2032.                                 if (++i >= ilim)
  2033.                                         break;
  2034.                                 eps *= 10.;
  2035.                                 d *= 10.;
  2036.                                 }
  2037.                         }
  2038.                 else {
  2039. #endif
  2040.                         /* Generate ilim digits, then fix them up. */
  2041.                         eps *= tens[ilim-1];
  2042.                         for(i = 1;; i++, d *= 10.) {
  2043.                                 L = (long)d;
  2044.                                 d -= L;
  2045.                                 *s++ = '0' + (int)L;
  2046.                                 if (i == ilim) {
  2047.                                         if (d > 0.5 + eps)
  2048.                                                 goto bump_up;
  2049.                                         else if (d < 0.5 - eps) {
  2050.                                                 while(*--s == '0');
  2051.                                                 s++;
  2052.                                                 goto ret1;
  2053.                                                 }
  2054.                                         break;
  2055.                                         }
  2056.                                 }
  2057. #ifndef No_leftright
  2058.                         }
  2059. #endif
  2060.  fast_failed:
  2061.                 s = s0;
  2062.                 d = d2;
  2063.                 k = k0;
  2064.                 ilim = ilim0;
  2065.                 }
  2066.  
  2067.         /* Do we have a "small" integer? */
  2068.  
  2069.         if (be >= 0 && k <= Int_max) {
  2070.                 /* Yes. */
  2071.                 ds = tens[k];
  2072.                 if (ndigits < 0 && ilim <= 0) {
  2073.                         if (ilim < 0 || d <= 5*ds)
  2074.                                 goto no_digits;
  2075.                         goto one_digit;
  2076.                         }
  2077.                 for(i = 1;; i++) {
  2078.                         L = (long)(d / ds);
  2079.                         d -= L*ds;
  2080. #ifdef Check_FLT_ROUNDS
  2081.                         /* If FLT_ROUNDS == 2, L will usually be high by 1 */
  2082.                         if (d < 0) {
  2083.                                 L--;
  2084.                                 d += ds;
  2085.                                 }
  2086. #endif
  2087.                         *s++ = '0' + (int)L;
  2088.                         if (i == ilim) {
  2089.                                 d += d;
  2090.                                 if (d > ds || d == ds && L & 1) {
  2091.  bump_up:
  2092.                                         while(*--s == '9')
  2093.                                                 if (s == s0) {
  2094.                                                         k++;
  2095.                                                         *s = '0';
  2096.                                                         break;
  2097.                                                         }
  2098.                                         ++*s++;
  2099.                                         }
  2100.                                 break;
  2101.                                 }
  2102.                         if (!(d *= 10.))
  2103.                                 break;
  2104.                         }
  2105.                 goto ret1;
  2106.                 }
  2107.  
  2108.         m2 = b2;
  2109.         m5 = b5;
  2110.         if (leftright) {
  2111.                 if (mode < 2) {
  2112.                         i =
  2113. #ifndef Sudden_Underflow
  2114.                                 denorm ? be + (Bias + (P-1) - 1 + 1) :
  2115. #endif
  2116. #ifdef IBM
  2117.                                 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
  2118. #else
  2119.                                 1 + P - bbits;
  2120. #endif
  2121.                         }
  2122.                 else {
  2123.                         j = ilim - 1;
  2124.                         if (m5 >= j)
  2125.                                 m5 -= j;
  2126.                         else {
  2127.                                 s5 += j -= m5;
  2128.                                 b5 += j;
  2129.                                 m5 = 0;
  2130.                                 }
  2131.                         if ((i = ilim) < 0) {
  2132.                                 m2 -= i;
  2133.                                 i = 0;
  2134.                                 }
  2135.                         }
  2136.                 b2 += i;
  2137.                 s2 += i;
  2138.                 mhi = i2b(Binit(&_mhi), 1);
  2139.                 }
  2140.         if (m2 > 0 && s2 > 0) {
  2141.                 i = m2 < s2 ? m2 : s2;
  2142.                 b2 -= i;
  2143.                 m2 -= i;
  2144.                 s2 -= i;
  2145.                 }
  2146.         if (b5 > 0) {
  2147.                 if (leftright) {
  2148.                         if (m5 > 0) {
  2149.                 Bigint *b_tmp;
  2150.                                 mhi = pow5mult(mhi, m5);
  2151.                                 b_tmp = mult(b_avail, mhi, b);
  2152.                                 b_avail = b;
  2153.                                 b = b_tmp;
  2154.                                 }
  2155.                         if (j = b5 - m5)
  2156.                                 b = pow5mult(b, j);
  2157.                         }
  2158.                 else
  2159.                         b = pow5mult(b, b5);
  2160.                 }
  2161.         S = i2b(S, 1);
  2162.         if (s5 > 0)
  2163.                 S = pow5mult(S, s5);
  2164.  
  2165.         /* Check for special case that d is a normalized power of 2. */
  2166.  
  2167.         if (mode < 2) {
  2168.                 if (!word1(d) && !(word0(d) & Bndry_mask)
  2169. #ifndef Sudden_Underflow
  2170.                  && word0(d) & Exp_mask
  2171. #endif
  2172.                                 ) {
  2173.                         /* The special case */
  2174.                         b2 += Log2P;
  2175.                         s2 += Log2P;
  2176.                         spec_case = 1;
  2177.                         }
  2178.                 else
  2179.                         spec_case = 0;
  2180.                 }
  2181.  
  2182.         /* Arrange for convenient computation of quotients:
  2183.          * shift left if necessary so divisor has 4 leading 0 bits.
  2184.          *
  2185.          * Perhaps we should just compute leading 28 bits of S once
  2186.          * and for all and pass them and a shift to quorem, so it
  2187.          * can do shifts and ors to compute the numerator for q.
  2188.          */
  2189.         if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f)
  2190.                 i = 32 - i;
  2191.         if (i > 4) {
  2192.                 i -= 4;
  2193.                 b2 += i;
  2194.                 m2 += i;
  2195.                 s2 += i;
  2196.                 }
  2197.         else if (i < 4) {
  2198.                 i += 28;
  2199.                 b2 += i;
  2200.                 m2 += i;
  2201.                 s2 += i;
  2202.                 }
  2203.         if (b2 > 0)
  2204.                 b = lshift(b, b2);
  2205.         if (s2 > 0)
  2206.                 S = lshift(S, s2);
  2207.         if (k_check) {
  2208.                 if (cmp(b,S) < 0) {
  2209.                         k--;
  2210.                         b = multadd(b, 10, 0);  /* we botched the k estimate */
  2211.                         if (leftright)
  2212.                                 mhi = multadd(mhi, 10, 0);
  2213.                         ilim = ilim1;
  2214.                         }
  2215.                 }
  2216.         if (ilim <= 0 && mode > 2) {
  2217.                 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
  2218.                         /* no digits, fcvt style */
  2219.  no_digits:
  2220.                         k = -1 - ndigits;
  2221.                         goto ret;
  2222.                         }
  2223.  one_digit:
  2224.                 *s++ = '1';
  2225.                 k++;
  2226.                 goto ret;
  2227.                 }
  2228.         if (leftright) {
  2229.                 if (m2 > 0)
  2230.                         mhi = lshift(mhi, m2);
  2231.  
  2232.                 /* Compute mlo -- check for special case
  2233.                  * that d is a normalized power of 2.
  2234.                  */
  2235.  
  2236.                 if (spec_case) {
  2237.             mlo = Brealloc(Binit(&_mlo), mhi->k);
  2238.                         Bcopy(mlo, mhi);
  2239.                         mhi = lshift(mhi, Log2P);
  2240.                         }
  2241.         else
  2242.             mlo = mhi;
  2243.  
  2244.                 for(i = 1;;i++) {
  2245.                         dig = quorem(b,S) + '0';
  2246.                         /* Do we yet have the shortest decimal string
  2247.                          * that will round to d?
  2248.                          */
  2249.                         j = cmp(b, mlo);
  2250.                         b_avail = diff(b_avail, S, mhi); /* b_avail = S - mi */
  2251.                         j1 = b_avail->sign ? 1 : cmp(b, b_avail);
  2252. #ifndef ROUND_BIASED
  2253.                         if (j1 == 0 && !mode && !(word1(d) & 1)) {
  2254.                                 if (dig == '9')
  2255.                                         goto round_9_up;
  2256.                                 if (j > 0)
  2257.                                         dig++;
  2258.                                 *s++ = dig;
  2259.                                 goto ret;
  2260.                                 }
  2261. #endif
  2262.                         if (j < 0 || j == 0 && !mode
  2263. #ifndef ROUND_BIASED
  2264.                                                         && !(word1(d) & 1)
  2265. #endif
  2266.                                         ) {
  2267.                                 if (j1 > 0) {
  2268.                                         b = lshift(b, 1);
  2269.                                         j1 = cmp(b, S);
  2270.                                         if ((j1 > 0 || j1 == 0 && dig & 1)
  2271.                                         && dig++ == '9')
  2272.                                                 goto round_9_up;
  2273.                                         }
  2274.                                 *s++ = dig;
  2275.                                 goto ret;
  2276.                                 }
  2277.                         if (j1 > 0) {
  2278.                                 if (dig == '9') { /* possible if i == 1 */
  2279.  round_9_up:
  2280.                                         *s++ = '9';
  2281.                                         goto roundoff;
  2282.                                         }
  2283.                                 *s++ = dig + 1;
  2284.                                 goto ret;
  2285.                                 }
  2286.                         *s++ = dig;
  2287.                         if (i == ilim)
  2288.                                 break;
  2289.                         b = multadd(b, 10, 0);
  2290.                         if (mlo == mhi)
  2291.                                 mlo = mhi = multadd(mhi, 10, 0);
  2292.                         else {
  2293.                                 mlo = multadd(mlo, 10, 0);
  2294.                                 mhi = multadd(mhi, 10, 0);
  2295.                                 }
  2296.                         }
  2297.                 }
  2298.         else
  2299.                 for(i = 1;; i++) {
  2300.                         *s++ = dig = quorem(b,S) + '0';
  2301.                         if (i >= ilim)
  2302.                                 break;
  2303.                         b = multadd(b, 10, 0);
  2304.                         }
  2305.  
  2306.         /* Round off last digit */
  2307.  
  2308.         b = lshift(b, 1);
  2309.         j = cmp(b, S);
  2310.         if (j > 0 || j == 0 && dig & 1) {
  2311.  roundoff:
  2312.                 while(*--s == '9')
  2313.                         if (s == s0) {
  2314.                                 k++;
  2315.                                 *s++ = '1';
  2316.                                 goto ret;
  2317.                                 }
  2318.                 ++*s++;
  2319.                 }
  2320.         else {
  2321.                 while(*--s == '0');
  2322.                 s++;
  2323.                 }
  2324.  ret:
  2325.     Bfree(b_avail);
  2326.         Bfree(S);
  2327.         if (mhi) {
  2328.                 if (mlo && mlo != mhi)
  2329.                         Bfree(mlo);
  2330.                 Bfree(mhi);
  2331.                 }
  2332.  ret1:
  2333.         Bfree(b);
  2334.         *s = 0;
  2335.         *decpt = k + 1;
  2336.         if (rve)
  2337.                 *rve = s;
  2338.         return s0;
  2339.         }
  2340. #endif /* USE_DTOA */
  2341.